setwd("") #Set working directory

# Read packages
library(haven)
library(list)
library(ggplot2)

# Read data
list_data <- read_dta("alt_estimators_data.dta")
list_data <- as.data.frame(list_data)

# Subset data
dir <- subset(list_data, dir_data == 1)

l_inc_list <- subset(list_data, h_inc1==1)
m_inc_list <- subset(list_data, h_inc2==1)
h_inc_list <- subset(list_data, h_inc3==1)

l_inc_dir <- subset(list_data, h_inc1==1 & dir_data == 1)
m_inc_dir <- subset(list_data, h_inc2==1 & dir_data == 1)
h_inc_dir <- subset(list_data, h_inc3==1 & dir_data == 1)

left_list <- subset(list_data, ideo_group1==1)
cent_list <- subset(list_data, ideo_group2==1)
right_list <- subset(list_data, ideo_group3==1)

left_dir <- subset(list_data, ideo_group1==1 & dir_data == 1)
cent_dir <- subset(list_data, ideo_group2==1 & dir_data == 1)
right_dir <- subset(list_data, ideo_group3==1 & dir_data == 1)

# Figure K.1
list_overall_ml <- ictreg(no_items ~ 1, data = list_data, J=4, method = "ml",overdispersed = TRUE)

dir_overall_ml <- glm(prog_dir ~ 1, data = dir, family = binomial("logit"))

sdb_all_ml <- predict(list_overall_ml, direct.glm = dir_overall_ml, se.fit = TRUE, level = 0.90)

sdb_all_ml_ <- data.frame(predicted_data = sdb_all_ml[1], z = 1:3, nam = c("List question", "Direct question", "Difference"))

ggplot(data = sdb_all_ml_, aes(x = factor(nam, level = c('List question', 'Direct question', 'Difference')), y = predicted_data.fit.fit)) +
  theme(
    plot.background = element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_rect(color="white",
                                fill=NA),
    panel.background = element_rect(fill = "white", colour = "grey50"),
    axis.line = element_line(colour = "black")) +
  geom_point() +
  geom_errorbar(aes(ymin = predicted_data.fit.lwr , ymax = predicted_data.fit.upr),width=0.1,
                size=0.3) +
  theme(plot.title = element_text(hjust = 0.5)) +
  ylab("Proportion opposing progressive taxation") +
  theme(legend.title = element_blank()) +
  theme(legend.position = "none") +
  geom_hline(yintercept=0,linetype=2) +
  theme(axis.text.x = element_text(color="black"),
        axis.text.y = element_text(color="black"),
        axis.title.x=element_blank(),
        axis.title.y = element_text(size = 8)) +
  theme(axis.ticks = element_line(colour = 'black')) +
  ylim(-0.5, 0.75)

# Figure K.2(a)
list_l_inc_ml <- ictreg(no_items ~ 1, data = l_inc_list, J=4, method = "ml",overdispersed = TRUE)

dir_l_inc_ml <- glm(prog_dir ~ 1, data = l_inc_dir, family = binomial("logit"))

sdb_l_inc_ml <- predict(list_l_inc_ml, direct.glm = dir_l_inc_ml, se.fit = TRUE, level = 0.90)

sdb_l_inc_ml_ <- data.frame(predicted_data = sdb_l_inc_ml[1], z = 1:3, nam = c("List question", "Direct question", "Difference"))

ggplot(data = sdb_l_inc_ml_, aes(x = factor(nam, level = c('List question', 'Direct question', 'Difference')), y = predicted_data.fit.fit)) +
  theme(
    plot.background = element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_rect(color="white",
                                fill=NA),
    panel.background = element_rect(fill = "white", colour = "grey50"),
    axis.line = element_line(colour = "black")) +
  geom_point() +
  geom_errorbar(aes(ymin = predicted_data.fit.lwr , ymax = predicted_data.fit.upr),width=0.1,
                size=0.3) +
  theme(plot.title = element_text(hjust = 0.5)) +
  ylab("Proportion opposing progressive taxation") +
  theme(legend.title = element_blank()) +
  theme(legend.position = "none") +
  geom_hline(yintercept=0,linetype=2) +
  theme(axis.text.x = element_text(color="black"),
        axis.text.y = element_text(color="black"),
        axis.title.x=element_blank(),
        axis.title.y = element_text(size = 8)) +
  theme(axis.ticks = element_line(colour = 'black')) +
  ylim(-0.5, 0.75)

# Figure K.2(b)
list_m_inc_ml <- ictreg(no_items ~ 1, data = m_inc_list, J=4, method = "ml",overdispersed = TRUE)

dir_m_inc_ml <- glm(prog_dir ~ 1, data = m_inc_dir, family = binomial("logit"))

sdb_m_inc_ml <- predict(list_m_inc_ml, direct.glm = dir_m_inc_ml, se.fit = TRUE, level = 0.90)

sdb_m_inc_ml_ <- data.frame(predicted_data = sdb_m_inc_ml[1], z = 1:3, nam = c("List question", "Direct question", "Difference"))

ggplot(data = sdb_m_inc_ml_, aes(x = factor(nam, level = c('List question', 'Direct question', 'Difference')), y = predicted_data.fit.fit)) +
  theme(
    plot.background = element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_rect(color="white",
                                fill=NA),
    panel.background = element_rect(fill = "white", colour = "grey50"),
    axis.line = element_line(colour = "black")) +
  geom_point() +
  geom_errorbar(aes(ymin = predicted_data.fit.lwr , ymax = predicted_data.fit.upr),width=0.1,
                size=0.3) +
  theme(plot.title = element_text(hjust = 0.5)) +
  ylab("Proportion opposing progressive taxation") +
  theme(legend.title = element_blank()) +
  theme(legend.position = "none") +
  geom_hline(yintercept=0,linetype=2) +
  theme(axis.text.x = element_text(color="black"),
        axis.text.y = element_text(color="black"),
        axis.title.x=element_blank(),
        axis.title.y = element_text(size = 8)) +
  theme(axis.ticks = element_line(colour = 'black')) +
  ylim(-0.5, 0.75)

# Figure K.2(c)
list_h_inc_ml <- ictreg(no_items ~ 1, data = h_inc_list, J=4, method = "ml",overdispersed = TRUE)

dir_h_inc_ml <- glm(prog_dir ~ 1, data = h_inc_dir, family = binomial("logit"))

sdb_h_inc_ml <- predict(list_h_inc_ml, direct.glm = dir_h_inc_ml, se.fit = TRUE, level = 0.90)

sdb_h_inc_ml_ <- data.frame(predicted_data = sdb_h_inc_ml[1], z = 1:3, nam = c("List question", "Direct question", "Difference"))

ggplot(data = sdb_h_inc_ml_, aes(x = factor(nam, level = c('List question', 'Direct question', 'Difference')), y = predicted_data.fit.fit)) +
  theme(
    plot.background = element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_rect(color="white",
                                fill=NA),
    panel.background = element_rect(fill = "white", colour = "grey50"),
    axis.line = element_line(colour = "black")) +
  geom_point() +
  geom_errorbar(aes(ymin = predicted_data.fit.lwr , ymax = predicted_data.fit.upr),width=0.1,
                size=0.3) +
  theme(plot.title = element_text(hjust = 0.5)) +
  ylab("Proportion opposing progressive taxation") +
  theme(legend.title = element_blank()) +
  theme(legend.position = "none") +
  geom_hline(yintercept=0,linetype=2) +
  theme(axis.text.x = element_text(color="black"),
        axis.text.y = element_text(color="black"),
        axis.title.x=element_blank(),
        axis.title.y = element_text(size = 8)) +
  theme(axis.ticks = element_line(colour = 'black')) +
  ylim(-0.5, 0.75)

# Figure K.3(a)
list_left_ml <- ictreg(no_items ~ 1, data = left_list, J=4, method = "ml",overdispersed = TRUE)

dir_left_ml <- glm(prog_dir ~ 1, data = left_dir, family = binomial("logit"))

sdb_left_ml <- predict(list_left_ml, direct.glm = dir_left_ml, se.fit = TRUE, level = 0.90)

sdb_left_ml_ <- data.frame(predicted_data = sdb_left_ml[1], z = 1:3, nam = c("List question", "Direct question", "Difference"))

ggplot(data = sdb_left_ml_, aes(x = factor(nam, level = c('List question', 'Direct question', 'Difference')), y = predicted_data.fit.fit)) +
  theme(
    plot.background = element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_rect(color="white",
                                fill=NA),
    panel.background = element_rect(fill = "white", colour = "grey50"),
    axis.line = element_line(colour = "black")) +
  geom_point() +
  geom_errorbar(aes(ymin = predicted_data.fit.lwr , ymax = predicted_data.fit.upr),width=0.1,
                size=0.3) +
  theme(plot.title = element_text(hjust = 0.5)) +
  ylab("Proportion opposing progressive taxation") +
  theme(legend.title = element_blank()) +
  theme(legend.position = "none") +
  geom_hline(yintercept=0,linetype=2) +
  theme(axis.text.x = element_text(color="black"),
        axis.text.y = element_text(color="black"),
        axis.title.x=element_blank(),
        axis.title.y = element_text(size = 8)) +
  theme(axis.ticks = element_line(colour = 'black')) +
  ylim(-0.5, 0.75)

# Figure K.3(b)
list_cent_ml <- ictreg(no_items ~ 1, data = cent_list, J=4, method = "ml",overdispersed = TRUE)

dir_cent_ml <- glm(prog_dir ~ 1, data = cent_dir, family = binomial("logit"))

sdb_cent_ml <- predict(list_cent_ml, direct.glm = dir_cent_ml, se.fit = TRUE, level = 0.90)

sdb_cent_ml_ <- data.frame(predicted_data = sdb_cent_ml[1], z = 1:3, nam = c("List question", "Direct question", "Difference"))

ggplot(data = sdb_cent_ml_, aes(x = factor(nam, level = c('List question', 'Direct question', 'Difference')), y = predicted_data.fit.fit)) +
  theme(
    plot.background = element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_rect(color="white",
                                fill=NA),
    panel.background = element_rect(fill = "white", colour = "grey50"),
    axis.line = element_line(colour = "black")) +
  geom_point() +
  geom_errorbar(aes(ymin = predicted_data.fit.lwr , ymax = predicted_data.fit.upr),width=0.1,
                size=0.3) +
  theme(plot.title = element_text(hjust = 0.5)) +
  ylab("Proportion opposing progressive taxation") +
  theme(legend.title = element_blank()) +
  theme(legend.position = "none") +
  geom_hline(yintercept=0,linetype=2) +
  theme(axis.text.x = element_text(color="black"),
        axis.text.y = element_text(color="black"),
        axis.title.x=element_blank(),
        axis.title.y = element_text(size = 8)) +
  theme(axis.ticks = element_line(colour = 'black')) +
  ylim(-0.5, 0.75)

# Figure K.3(c)
list_right_ml <- ictreg(no_items ~ 1, data = right_list, J=4, method = "ml",overdispersed = TRUE)

dir_right_ml <- glm(prog_dir ~ 1, data = right_dir, family = binomial("logit"))

sdb_right_ml <- predict(list_right_ml, direct.glm = dir_right_ml, se.fit = TRUE, level = 0.90)

sdb_right_ml_ <- data.frame(predicted_data = sdb_right_ml[1], z = 1:3, nam = c("List question", "Direct question", "Difference"))

ggplot(data = sdb_right_ml_, aes(x = factor(nam, level = c('List question', 'Direct question', 'Difference')), y = predicted_data.fit.fit)) +
  theme(
    plot.background = element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_rect(color="white",
                                fill=NA),
    panel.background = element_rect(fill = "white", colour = "grey50"),
    axis.line = element_line(colour = "black")) +
  geom_point() +
  geom_errorbar(aes(ymin = predicted_data.fit.lwr , ymax = predicted_data.fit.upr),width=0.1,
                size=0.3) +
  theme(plot.title = element_text(hjust = 0.5)) +
  ylab("Proportion opposing progressive taxation") +
  theme(legend.title = element_blank()) +
  theme(legend.position = "none") +
  geom_hline(yintercept=0,linetype=2) +
  theme(axis.text.x = element_text(color="black"),
        axis.text.y = element_text(color="black"),
        axis.title.x=element_blank(),
        axis.title.y = element_text(size = 8)) +
  theme(axis.ticks = element_line(colour = 'black')) +
  ylim(-0.5, 0.75)

# Figure K.4
list_overall_nls <- ictreg(no_items ~ 1, data = list_data, J=4, method = "nls")

dir_overall_nls <- glm(prog_dir ~ 1, data = dir, family = binomial("logit"))

sdb_all_nls <- predict(list_overall_nls, direct.glm = dir_overall_nls, se.fit = TRUE, level = 0.90)

sdb_all_nls_ <- data.frame(predicted_data = sdb_all_nls[1], z = 1:3, nam = c("List question", "Direct question", "Difference"))

ggplot(data = sdb_all_nls_, aes(x = factor(nam, level = c('List question', 'Direct question', 'Difference')), y = predicted_data.fit.fit)) +
  theme(
    plot.background = element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_rect(color="white",
                                fill=NA),
    panel.background = element_rect(fill = "white", colour = "grey50"),
    axis.line = element_line(colour = "black")) +
  geom_point() +
  geom_errorbar(aes(ymin = predicted_data.fit.lwr , ymax = predicted_data.fit.upr),width=0.1,
                size=0.3) +
  theme(plot.title = element_text(hjust = 0.5)) +
  ylab("Proportion opposing progressive taxation") +
  theme(legend.title = element_blank()) +
  theme(legend.position = "none") +
  geom_hline(yintercept=0,linetype=2) +
  theme(axis.text.x = element_text(color="black"),
        axis.text.y = element_text(color="black"),
        axis.title.x=element_blank(),
        axis.title.y = element_text(size = 8)) +
  theme(axis.ticks = element_line(colour = 'black')) +
  ylim(-0.5, 0.75)

# Figure K.5(a)
list_l_inc_nls <- ictreg(no_items ~ 1, data = l_inc_list, J=4, method = "nls")

dir_l_inc_nls <- glm(prog_dir ~ 1, data = l_inc_dir, family = binomial("logit"))

sdb_l_inc_nls <- predict(list_l_inc_nls, direct.glm = dir_l_inc_nls, se.fit = TRUE, level = 0.90)

sdb_l_inc_nls_ <- data.frame(predicted_data = sdb_l_inc_nls[1], z = 1:3, nam = c("List question", "Direct question", "Difference"))

ggplot(data = sdb_l_inc_nls_, aes(x = factor(nam, level = c('List question', 'Direct question', 'Difference')), y = predicted_data.fit.fit)) +
  theme(
    plot.background = element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_rect(color="white",
                                fill=NA),
    panel.background = element_rect(fill = "white", colour = "grey50"),
    axis.line = element_line(colour = "black")) +
  geom_point() +
  geom_errorbar(aes(ymin = predicted_data.fit.lwr , ymax = predicted_data.fit.upr),width=0.1,
                size=0.3) +
  theme(plot.title = element_text(hjust = 0.5)) +
  ylab("Proportion opposing progressive taxation") +
  theme(legend.title = element_blank()) +
  theme(legend.position = "none") +
  geom_hline(yintercept=0,linetype=2) +
  theme(axis.text.x = element_text(color="black"),
        axis.text.y = element_text(color="black"),
        axis.title.x=element_blank(),
        axis.title.y = element_text(size = 8)) +
  theme(axis.ticks = element_line(colour = 'black')) +
  ylim(-0.5, 0.75)

# Figure K.5(b)
list_m_inc_nls <- ictreg(no_items ~ 1, data = m_inc_list, J=4, method = "nls")

dir_m_inc_nls <- glm(prog_dir ~ 1, data = m_inc_dir, family = binomial("logit"))

sdb_m_inc_nls <- predict(list_m_inc_nls, direct.glm = dir_m_inc_nls, se.fit = TRUE, level = 0.90)

sdb_m_inc_nls_ <- data.frame(predicted_data = sdb_m_inc_nls[1], z = 1:3, nam = c("List question", "Direct question", "Difference"))

ggplot(data = sdb_m_inc_nls_, aes(x = factor(nam, level = c('List question', 'Direct question', 'Difference')), y = predicted_data.fit.fit)) +
  theme(
    plot.background = element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_rect(color="white",
                                fill=NA),
    panel.background = element_rect(fill = "white", colour = "grey50"),
    axis.line = element_line(colour = "black")) +
  geom_point() +
  geom_errorbar(aes(ymin = predicted_data.fit.lwr , ymax = predicted_data.fit.upr),width=0.1,
                size=0.3) +
  theme(plot.title = element_text(hjust = 0.5)) +
  ylab("Proportion opposing progressive taxation") +
  theme(legend.title = element_blank()) +
  theme(legend.position = "none") +
  geom_hline(yintercept=0,linetype=2) +
  theme(axis.text.x = element_text(color="black"),
        axis.text.y = element_text(color="black"),
        axis.title.x=element_blank(),
        axis.title.y = element_text(size = 8)) +
  theme(axis.ticks = element_line(colour = 'black')) +
  ylim(-0.5, 0.75)

# Figure K.5(c)
list_h_inc_nls <- ictreg(no_items ~ 1, data = h_inc_list, J=4, method = "nls")

dir_h_inc_nls <- glm(prog_dir ~ 1, data = h_inc_dir, family = binomial("logit"))

sdb_h_inc_nls <- predict(list_h_inc_nls, direct.glm = dir_h_inc_nls, se.fit = TRUE, level = 0.90)

sdb_h_inc_nls_ <- data.frame(predicted_data = sdb_h_inc_nls[1], z = 1:3, nam = c("List question", "Direct question", "Difference"))

ggplot(data = sdb_h_inc_nls_, aes(x = factor(nam, level = c('List question', 'Direct question', 'Difference')), y = predicted_data.fit.fit)) +
  theme(
    plot.background = element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_rect(color="white",
                                fill=NA),
    panel.background = element_rect(fill = "white", colour = "grey50"),
    axis.line = element_line(colour = "black")) +
  geom_point() +
  geom_errorbar(aes(ymin = predicted_data.fit.lwr , ymax = predicted_data.fit.upr),width=0.1,
                size=0.3) +
  theme(plot.title = element_text(hjust = 0.5)) +
  ylab("Proportion opposing progressive taxation") +
  theme(legend.title = element_blank()) +
  theme(legend.position = "none") +
  geom_hline(yintercept=0,linetype=2) +
  theme(axis.text.x = element_text(color="black"),
        axis.text.y = element_text(color="black"),
        axis.title.x=element_blank(),
        axis.title.y = element_text(size = 8)) +
  theme(axis.ticks = element_line(colour = 'black')) +
  ylim(-0.5, 0.75)

# Figure K.6(a)
list_left_nls <- ictreg(no_items ~ 1, data = left_list, J=4, method = "nls")

dir_left_nls <- glm(prog_dir ~ 1, data = left_dir, family = binomial("logit"))

sdb_left_nls <- predict(list_left_nls, direct.glm = dir_left_nls, se.fit = TRUE, level = 0.90)

sdb_left_nls_ <- data.frame(predicted_data = sdb_left_nls[1], z = 1:3, nam = c("List question", "Direct question", "Difference"))

ggplot(data = sdb_left_nls_, aes(x = factor(nam, level = c('List question', 'Direct question', 'Difference')), y = predicted_data.fit.fit)) +
  theme(
    plot.background = element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_rect(color="white",
                                fill=NA),
    panel.background = element_rect(fill = "white", colour = "grey50"),
    axis.line = element_line(colour = "black")) +
  geom_point() +
  geom_errorbar(aes(ymin = predicted_data.fit.lwr , ymax = predicted_data.fit.upr),width=0.1,
                size=0.3) +
  theme(plot.title = element_text(hjust = 0.5)) +
  ylab("Proportion opposing progressive taxation") +
  theme(legend.title = element_blank()) +
  theme(legend.position = "none") +
  geom_hline(yintercept=0,linetype=2) +
  theme(axis.text.x = element_text(color="black"),
        axis.text.y = element_text(color="black"),
        axis.title.x=element_blank(),
        axis.title.y = element_text(size = 8)) +
  theme(axis.ticks = element_line(colour = 'black')) +
  ylim(-0.5, 0.75)

# Figure K.6(b)
list_cent_nls <- ictreg(no_items ~ 1, data = cent_list, J=4, method = "nls")

dir_cent_nls <- glm(prog_dir ~ 1, data = cent_dir, family = binomial("logit"))

sdb_cent_nls <- predict(list_cent_nls, direct.glm = dir_cent_nls, se.fit = TRUE, level = 0.90)

sdb_cent_nls_ <- data.frame(predicted_data = sdb_cent_nls[1], z = 1:3, nam = c("List question", "Direct question", "Difference"))

ggplot(data = sdb_cent_nls_, aes(x = factor(nam, level = c('List question', 'Direct question', 'Difference')), y = predicted_data.fit.fit)) +
  theme(
    plot.background = element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_rect(color="white",
                                fill=NA),
    panel.background = element_rect(fill = "white", colour = "grey50"),
    axis.line = element_line(colour = "black")) +
  geom_point() +
  geom_errorbar(aes(ymin = predicted_data.fit.lwr , ymax = predicted_data.fit.upr),width=0.1,
                size=0.3) +
  theme(plot.title = element_text(hjust = 0.5)) +
  ylab("Proportion opposing progressive taxation") +
  theme(legend.title = element_blank()) +
  theme(legend.position = "none") +
  geom_hline(yintercept=0,linetype=2) +
  theme(axis.text.x = element_text(color="black"),
        axis.text.y = element_text(color="black"),
        axis.title.x=element_blank(),
        axis.title.y = element_text(size = 8)) +
  theme(axis.ticks = element_line(colour = 'black')) +
  ylim(-0.5, 0.75)

# Figure K.6(c)
list_right_nls <- ictreg(no_items ~ 1, data = right_list, J=4, method = "nls")

dir_right_nls <- glm(prog_dir ~ 1, data = right_dir, family = binomial("logit"))

sdb_right_nls <- predict(list_right_nls, direct.glm = dir_right_nls, se.fit = TRUE, level = 0.90)

sdb_right_nls_ <- data.frame(predicted_data = sdb_right_nls[1], z = 1:3, nam = c("List question", "Direct question", "Difference"))

ggplot(data = sdb_right_nls_, aes(x = factor(nam, level = c('List question', 'Direct question', 'Difference')), y = predicted_data.fit.fit)) +
  theme(
    plot.background = element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_rect(color="white",
                                fill=NA),
    panel.background = element_rect(fill = "white", colour = "grey50"),
    axis.line = element_line(colour = "black")) +
  geom_point() +
  geom_errorbar(aes(ymin = predicted_data.fit.lwr , ymax = predicted_data.fit.upr),width=0.1,
                size=0.3) +
  theme(plot.title = element_text(hjust = 0.5)) +
  ylab("Proportion opposing progressive taxation") +
  theme(legend.title = element_blank()) +
  theme(legend.position = "none") +
  geom_hline(yintercept=0,linetype=2) +
  theme(axis.text.x = element_text(color="black"),
        axis.text.y = element_text(color="black"),
        axis.title.x=element_blank(),
        axis.title.y = element_text(size = 8)) +
  theme(axis.ticks = element_line(colour = 'black')) +
  ylim(-0.5, 0.75)